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We evaluate the probability that a Boolean network returns to an attractor after perturbing 
h nodes. We find that the return probability as function of h can display a variety of different 
behaviours, which yields insights into the state-space structure. In addition to performing computer 
simulations, we derive analytical results for several types of Boolean networks, in particular for 
Random Boolean Networks. We also apply our method to networks that have been evolved for 
robustness to small perturbations, and to a biological example. 
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I. INTRODUCTION 

Boolean networks (BN) are used to model the dynam- 
ics of a wide variety of complex systems, ranging from 
neural networks and social systems to gene regu- 
lation networks [j], sometimes combined with evolution- 
ary processes Q ■ BN are composed of interacting nodes 
with binary states, typically and 1, coupled among each 
other. The state of each node evolves according to a func- 
tion of the states observed in a certain neighbourhood, 
similar to what is done when using cellular automata Q , 
but in contrast to cellular automata, BN have no regular 
lattice structure, and not all nodes are assigned the same 
update function. 

Usually, BN models are studied using deterministic dy- 
namical rules. After a transient time, these networks 
reach attractors, which are a sequence of periodically re- 
peated states. The number and length of such attractors 
is an important property investigated in BNs. Of equal 
importance are the sizes of the basins of attraction of 
these attractors, which are the sets of states leading to 
the attractors. Some networks may contain periodic se- 
quences in state space that are not "attractors" in the 
strict sense, because there are no states outside this se- 
quence that are attracted to it. However, in this paper 
we do not make this distinction, and we call all sequence 
periodically repeated sequences of states "attractor" . 

Real networks are often influenced by noise, since 
molecule concentrations may be small (e.g., in biologi- 
cal systems 0]) or behaviour may be unpredictable (e.g., 
in social systems). For this reason, it is important to 
investigate how robust the behaviour found under deter- 
ministic update rules is when noise is added. Examples 
for such studies are Ising-like models placed on a network 
topology 0, 3 ; BN with a probabilistic rule for choosing 
the update function at each time stepQ, models with 
stochastic update sequences [l^ and small stochastic de- 
lays in the update time [ll| . These stochastic models can 
lead to surprising new results. For instance, for stochas- 
tic update sequences it could be shown that the number 
of attractors (now defined as a recurrent set of states) 
grows like a power law [l^, [U as function of the net- 
work size, while it grows superpolynoniially for parallel 
update. 



In this paper, we investigate the effect of a perturba- 
tion on a BN that is updated in parallel and that is on 
an attractor. The quantity we evaluate is the probability 
that the network returns to the same attractor after the 
perturbation, as a function of the size of the perturba- 
tion, i.e., of the number of nodes the states of which are 
changed. This leads to a curve that is characteristic of 
the network and is strikingly different for different types 
of networks. 

The outline of this paper is as follows: First, we investi- 
gate two types of simple networks, namely networks con- 
sisting of independent nodes and networks consisting of a 
single loop of nodes. The results then help to understand 
the behaviour of Random Boolean Networks under per- 
turbations, which are studied in Section [IIII dealing with 
frozen, critical and chaotic networks. In Section llVl we 
then investigate a few specific networks, which are not 
random networks, and we find that their characteristic 
curves are very different from those of random networks, 
reflecting for instance the higher robustness to perturba- 
tions. Finally, we summarize and discuss our findings in 
Section IVl 



II. SIMPLE NETWORKS 

A. Independent nodes 

Let us first consider a system of N independent nodes. 
Then the response of each node to the perturbation is 
independent of the response of the other nodes. The 
state of each of the N nodes at time t -I- 1 is determined 
by its state at time t. There are 4 possibilities to assign to 
such a node an update function: (i) the state of the node 
is irrespective of the state at the previous time step; 
(ii) the state of the node is 1 irrespective of the state at 
the previous time step; (iii) the state of the node at time 
t 4- 1 is identical to the state at time t; (iv) the state of 
the node at time < -I- 1 is the opposite of the state at time 
t. This means that the node alternates between and 
1. The first two update functions are constant functions, 
the third is the "copy" function, the fourth the "invert" 
function. 

We generate a network by assigning update functions 
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to the TV nodes. We then initiahze the network in a ran- 
domly chosen state and wait until it reaches an attractor. 
In the simple networks considered here, an attractor is 
reached after one time step. Unless no "invert" function 
is chosen, the period of the attractor is 2. When a node is 
perturbed, its response depends on its update function: 
If the update function is constant, the node returns in 
the next time step to its previous value. If the update 
function is "copy" , the node remains in its new state and 
does not return to its previous state. If the update func- 
tion is "invert" , the node continues to oscillate between 
1 and 0, however with a phase shift of one time step. 

When we perturb h nodes, the network returns to the 
attractor only if all perturbed nodes have a constant func- 
tion, or if all nodes with an "invert" function are per- 
turbed but none with a "copy" function. Let us focus 
on the case that N is large and that each of the four 
update functions is assigned to one quarter of the nodes. 
The probability that the network returns to its previous 
attractor after perturbing h nodes is then 

(N/2\ 

Pret{h) = ^, (1) 
\h) 

which can be approximated by 

Pretih) ^ (2) 



The most important result of this subsection is that 
for small perturbations the return probability is simply 
given by an exponential dependence 

Pretih) ^ Pretil)'' ■ 

This result will be generally true for small h whenever 
perturbations at different nodes decay independently 
from each other. For this reason, we will see an expo- 
nential decrease in the examples below when N is large 
and h is small. 



B. Simple loops 

Now, let us consider a simple loop of N nodes. 




for h/N <^ 1. For h > N/2, the return probability is 
0, since it is no longer possible that all perturbed nodes 
have constant hmctions. For N/A < h < 3N/4, it is in 
principle possible that all nodes with "invert" functions 
but none with "copy" functions are perturbed, however, 
the probability that this occurs is so small that we neglect 
it here. 

For general networks, the return probability has to be 
evaluated by averaging over different initial states, how- 
ever, in the special case considered here the return prob- 
ability is the same for each initial state. 
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FIG. 1: Pretih) (in percent) of a set of independent nodes as 
given by Eq. ((TJ, the dashed line is obtained using Eq. © 



FIG. 2: A loop of seven nodes 

Each node is connected to its predecessor on the loop, 
and the update function of each node is either "copy" 
or "invert" . A constant function at one node will freeze 
the whole loop, in which case Pret = 1 for all h, and we 
therefore focus on the more interesting case of no con- 
stant function in the loop. A loop with n inversions can 
be mapped bijectively onto one with (n — 2) inversions 
by replacing two "invert" with two "copy" functions and 
by inverting the values of all nodes between these two 
couplings. It is therefore sufRcient to distinguish loops 
with an even or an odd number of inversions, and we call 
them "even" and "odd" loops respectively. When dis- 
cussing even loops, we consider loops with only "copy" 
functions. To odd loops we assign one "invert" function 
and — 1 "copy" functions. An even loop with a prime 
number of nodes returns to its initial state after A^ time 
steps. If N is not prime, shorter periods exist. Further- 
more, if all nodes have the same state, the loop is on a 
fixed point. An odd loop with a prime number of nodes 
returns to its initial state after 2A^ time steps. An odd 
loop has no fixed points. The shortest attractor has pe- 
riod 2, with alternating I's and O's. 

Figure [3] shows Pret for an even loop with A^ = 7, as 
obtained from computer simulations. 
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FIG. 3; Pret (ft) (in percent) for an even loop with N — 7 



The return probability is for all odd values of ft., 
and it has the same nonzero value for all even h larger 
than 0. This happens for all even loops with N being a 
prime number, as we show in the following. If h nodes 
are perturbed, the loop returns to the same attractor 
only if the perturbation leads it to another state of this 
attractor. All states of the same attractor have the same 
sequence of Os and Is, but rotated around the loop by 
some number of steps. Only by inverting an identical 
number of Os and Is can we stay on the same attractor, 
and this explains why Pret is nonzero only for even h. 
Now, by fixing the nodes that are to be perturbed and by 
fixing the number m of steps by which the perturbation 
of these nodes shall rotate the attractor, we uniquely fix 
the attractor (apart from an inversion of all nodes) . This 
attractor is found by fixing the state of one node to 1 (or 
0) and by requiring that the node m steps ahead has the 
same state if the first node is not a perturbed node, and 
the opposite state if the first node is a perturbed node. 
Then wc determine in the same way the state of the node 
m steps further, and so on, until the state of all nodes is 
fixed. Since A'^ is a prime number, we do not return to 
the first node before all other nodes have been visited. 
The probability to be still on the same attractor after 
perturbing h nodes is for even h therefore 



_ 2(7V-1) 



which is the number of configurations that are rotated by 
m = 1, . . . iV — 1 steps under the perturbation, divided by 
the total number of configurations of the loop. If N is not 
a prime number, there are attractors of different length, 
and the function PretW depends in a more complicated 
way on N . 
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FIG. 4: Pret (ft) (in percent) of an odd loop N = 7 

For an odd loop with prime N, Pret is 1 for /i = and 
h = N, and has otherwise the value Pret = ^{N ~ l)/2^. 
This result is obtained in a similar way as for even loops: 
The loop remains on the same attractor after perturb- 
ing h nodes, if the state generated by the perturbation is 
among the 27V — 2 different states assumed by the loop 
during the next 2A^ — 1 time steps, but not after N steps 
(where the state of all nodes is simply inverted, which can 
only be achieved for h = N). If we fix the nodes to be 
perturbed and the number m of time steps between the 
unperturbed and the perturbed state, we can uniquely 
identify the attractor by fixing the state of one node and 
then stepping around the loop in steps of size m, assign- 
ing to each node the same state as to the previous node 
if the previous node is not among the perturbed nodes 
and if the "invert" function is not between the two nodes. 
Otherwise, the inverted state is assigned. 



C. Collections of several loops 

If the network consists of several independent loops, 
the periods of the attractors are the least common mul- 
tiples of the periods of the loops. Obviously, a collection 
of loops can remain on the same attractor after a per- 
turbation only if the perturbation produces a state that 
corresponds to a future state of the loop, but not a state 
that is not part of the state sequence on this loop. If all 
loop lengths are different from each other and have no 
common divisor, every combination of future loop states 
belongs to the same attractor, and the return probabil- 
ity is identical to the probability that the perturbation 
generates only such states. Otherwise, if there are two or 
more loops with periods that have a common divisor, the 
perturbation must advance all these loops by a multiple 
of this divisor if the system shall remain on the same at- 
tractor. Since the analytical expressions for Pret become 
complicated and do not provide special insight, we omit 
them here. We only note that the function Pret{h) de- 
pends strongly on the lengths and number of loops, and 
on the common divisors of their lengths. 



4 



III. RANDOM BOOLEAN NETWORKS 

A random Boolean network (RBN) is constructed by 
choosing for each node at random k nodes from which it 
receives its input, and an update function that assigns to 
each of the 2'^ states of the k input nodes an output 1 or 0. 
The update function of each node is chosen at random 
among all 2^ possible update functions [H, [l^. All 
nodes are updated in parallel. Depending on the value of 
k and the probabilities assigned to the different update 
functions, the dynamics of the network is either in the 
frozen or in the chaotic phase. At the boundary between 
the two are critical networks. If all update functions are 
assigned the same probability, RBNs with fc = 1 are in 
the frozen phase, networks with k = 2 are critical, and 
networks with k > 2 are chaotic. In the following, we 
consider perturbations of these three types of networks. 

A. RBNs in the frozen phase 

In the frozen phase, all nodes apart from a small num- 
ber (that remains finite in the limit of infinite system 
size) assume a constant value after a transient time. If in 
the stationary state the value of one node is changed, this 
perturbation propagates during one time step on average 
to less than one other node. If all nodes that become 
frozen after some time arc removed from the network, 
there remain the nonfrozcn nodes. In the limit N oo, 
it can be shown that these nonfrozen nodes are connected 
to simple loops (the "relevant loops") with trees rooted 
in the loops [l^ . The nodes in the trees are slaved to the 
dynamics on the loops, and the perturbation of a node 
on a tree does not induce a change of the attractor. The 
perturbation of a frozen node can affect a couple of other 
nodes to which the perturbation may propagate. If these 
other nodes are also frozen or part of a nonfrozen tree, 
they will soon return to the behaviour they showed before 
the perturbation. However, if the perturbation affects a 
node sitting on a relevant loop, the attractor will usually 
be changed. If we denote with pr the probability that 
the perturbation of one node will affect a node on a rel- 
evant loop, we obtain Pret{h) — (1 — Pr)^ for small h. In 
the limit N oo, the probability pr must become pro- 
portional to 1/N, since the number of nodes on relevant 
loops remains finite in this limit. 

From these properties of frozen networks in the limit 
of large N, we can conclude that Pret{h) is close to 1 
for h <^ N, since the probability of perturbing one of 
the relevant nodes is very small. When h becomes of 
the order of Npr, relevant nodes are perturbed with a 
considerable probability, and the shape of the function 
Pret{h) depends on the properties of the relevant loops. 
The size and number of relevant loops is in general dif- 
ferent in different networks. From our discussion in the 
previous section for collection of independent loops, we 
conclude that there is no self-averaging of P^etih) in the 
limit of large N . Also, the probability that the number 



of relevant nodes is approaches a nonzero constant in 
the limit N ^ oo, and therefore some networks always 
have a constant return probability Pret{h) = 1. 

Figure [5] shows the result of a computer simulation of 
three networks in the frozen phase. Half of the nodes 
were assigned the value k = 1 and half of the nodes the 
value k = 2. (This combination of k values was chosen 
because networks with fc = 1 are almost always com- 
pletely frozen. Even for the value (fc*) — 1.5, only few 
networks were not completely frozen.) One can see the 
inital exponential decay and the broad range of curve 
shapes when h/N is not small anymore. The behaviour 
of the curves when h approaches N depends on the prob- 
ability that the inversion of the state of all nodes leaves 
the network in the same basin of attraction. As wc can 
see in this and other figures in this paper, this probability 
varies widely between networks and is not rarely larger 
than Pret{N/2). 
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FIG. 5: Pret{h)(in percent) of three frozen networks = 
200, (fc> = 1.5 



B. RBNs in the chaotic phase 

In the chaotic phase, initially similar configurations di- 
verge exponentially. Attractors are usually long, and a 
non- vanishing proportion of all nodes keep changing their 
state even after long times. Networks with k = N are 
easiest to understand among the chaotic networks, since 
the successor of each network state is a randomly chosen 
other network state [lB|. In state space, we have there- 
fore a "network" with one randomly chosen "input" (i.e. 
successor) for each "node" (i.e., state). In such a k = N 
network, every perturbation leads the network anywhere 
in state space, and wc therefore expect Pret{h) to be a 
constant function, with the exception of the value 1 at 
/i = 0. The average value of this constant can be de- 
termined from what is known for k = N networks: The 
value of Pret on the plateau is given by the average proba- 
bility that a randomly chosen state in state space belongs 
to the same basin of attraction as the randomly chosen 
initial network state. We define the weight Wp of an at- 
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tractor p as the length of the attractor plus the number 
of basin states draming into that attractor, normalized 
by the size of the state space (2^), so that X]p ^ 1 
□3. We therefore have 



Pr. 



(4) 



In the limit N ^ oo, the average number of attractors 
with a weight between w and w + dw is given by [igI . [l7j 



9{w) 



2w\/l — w ' 



leading to 



g{w)w dw = — 
3 



(5) 



(6) 



This analytical result applies to the ensemble average of 
all networks with k = N . In a given network, the value 
of the plateau usually deviates from 2/3, since always 
a few large basins dominate the value of Pret, and the 
values w of these basins differ between networks even in 
the thermodynamic limit. 
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FIG. 6: Pret{h){m percent) of three chaotic networks with 
TV = 60,fc = 3 

Now, let us turn to chaotic networks with fixed k 
(smaller than N), for which no such analytical results 
are known. Figure [6] shows Pret as function of h for three 
networks with A; = 3 and = 60. We see a plateau for 
h > 2, which means that perturbing one or two nodes 
does not yet necessarily carry the network to a random 
state, but perturbing several nodes does. The reason for 
this is that a certain proportion of nodes in chaotic net- 
works with fixed k assume a constant value after some 
time, and perturbing only such nodes may not change 
the attractor. For larger h, it becomes very likely that 
relevant nodes are changed by a perturbation, and then 
the situation is similar as in fc = networks, which have 
no frozen core. 



Chaotic networks with fixed value of k and large N 
have been shown to share many properties with chaotic 
k ~ N networks. From [l^, it appears that chaotic net- 
works with k < N have the same set of basin weights as 
those with k ^ N, which would mean that the average 
height of the plateau should be the same for all chaotic 
networks in the thermodynamic limit. Wc tested this as- 
sumption by simulating smaller networks (of size 20) on 
the computer for fc = 3, 4 and 5. The average of the 
return probability for a perturbation of size N/2 = 10 
over approximately 10000 networks is Pret(lO) — 0.6764, 
0.675, 0.6735, 0.6686 and 0.6675 for k = 3,4,5,6,9 re- 
spectively. These values are impressively close to 2/3 
given the small system size. For larger fc, the value is 
closer to 2/3, because the number of nonfrozcn nodes 
and therefore the size of the (nontrivial part of the) state 
space increases with increasing k. We also simulated sys- 
tems with larger values of N, however, it is very hard to 
obtain good statistics for larger N because the size of the 
state space and the length of attractors increase exponen- 
tially with N. From the data we have it appears that the 
value of the plateau moves closer to 2/3 with increas- 
ing N. Figure [7] shows the probability distributions of 
Pret (10), from which the above-mentioned averages have 
been obtained. This distribution is very broad and does 
not appear to become narrower for larger N. There is no 
self-averaging of the plateau value of Pret ■ 

In summary, we find that chaotic networks with large 
N have a return probability Pret{h) that decays rapidly 
for small h and then reaches a plateau, the ensemble 
average of which lies at 2/3. 
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FIG. 7: Probablity distribution of the value Pret(lO) for ap- 
proximately 10000 networks of size A'' = 20, for different val- 
ues of k. 



C. Critical RBNs 

Critical networks are at the boundary between the 
frozen and chaotic phase, and neighbouring configura- 
tions diverge only algebraically with time. Just as frozen 
networks, they have a large frozen core. For large iV, the 
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number of nodes that do not become frozen after some 
time, is proportional to N'^/^ [l^. Most of them are ir- 
relevant for determining the attractors. When the frozen 
core is removed from the network, the remaining nodes 
are connected to relevant components, most of which are 
simple loops and all of which contain loops, and with 
trees rooted in these loops. The number of nodes in the 
relevant components is proportional to N^^'^ for large N . 
The mean number of nodes affected by the perturbation 
of a single node is ~ N'^/^ for large N . The probabil- 
ity that a perturbation of size h = 1 affects a relevant 
node is given by the probability ~ 7V~^/^ that a given 
node is affected by the perturbation, times the number of 
relevant nodes, Prei(l) — aN^^^-^. We therefore obtain 

Pretih) ~ (a7V-i/3)" (7) 

for h <^ iV^/^. This exponential decay can be seen in 
Figure [5] for small h. For h > iV^/^, the probabihty of 
perturbing a relevant node is not small anymore, and P^et 
reaches a plateau. A rough estimate of the dependence 
on A'' of the value of Pret on the plateau is obtained by the 
following reasoning: The number of attractors increases 
with N roughly as 2^ '' , and the basins of attraction 
have a size of the order 2^/2^'''^ Using Eq. (g]), this 

gives the estimate Pret ^ 2~^^^'^ on the plateau with 
some constant h. The height of the plateau decreases 
rapidly with increasing N because of the vast number of 
different attractors of comparable basin size. 



critical networks. For such small networks sizes, the crit- 
ical features of Pretih) derived for the thermodynamic 
limit cannot yet be seen. 



IV. SOME SPECIFIC NETWORKS 

In the ensemble of all random networks of a not too 
large size, there are some networks whose characteristic 
curves are very different from conventional ones, as shown 
in the next two figures. 
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FIG. 9: Pret{h){m percent) of a critical network with A'' = 
100, fc = 2 
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FIG. 8: Pret (ft) (in percent) of three critical networks A'" = 
100, fc = 2 

Taking these results together, we expect that in large 
critical networks the function Pretih) decays from 1 to a 
value close to when h increases from to N^^^ and then 
stays on this plateau as h increases further. The value of 
plateau decreases as an exponential function of —N^^"^. 
However, it is known that the scaling with N'^/^ and N^^^ 
of the nonfrozen and relevant nodes becomes clearly visi- 
ble only in huge networks (of the order 10^), and smaller 
networks may show broad distributions in the numbers 
of these nodes. Figure [S] shows three examples of Pret for 
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FIG. 10: Pret (ft) (in percent) of another critical network with 
iV = 300, k = 2 

The first curve is that of a very robust network, which 
additionally has the property of returning more easily to 
the original attractor when the perturbation is larger. 

The second curve is that of a network that has no 
plateau but a very extended decrease of Pret with h. A 
closer inspection of this network reveals that, although 
its parameter values classify it as critical, there is no 
frozen core but all nodes are relevant and part of a single 
complex relevant component. This means that the state 
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space is far from random, and that perturbations of dif- 
ferent sizes carry the network to different regions in state 
space. 

Let us compare the first curve to that of networks 
that have been evolved for robustness to small pertur- 
bations [13 , shown in Figure [TlJ These networks have 
been evolved by modifying the connections and functions 
and by selecting for a high probability of returning to 
the same attractor after perturbing one node. Interest- 
ingly, although only a large value of Preti^) has been im- 
posed in the system, Pret is very high for all perturbation 
sizes. Many networks obtained by the same procedure 
even show Pret{h) = 1 for all h. 
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FIG. 11: Pret{h.){in percent) of a network of 60 nodes evolved 
for robustness. 



The network used for producing Figure [12] is the core 
part of cell cycle regulation network of budding yeast, as 
represented in [20| . An important property of this net- 
work is that it has a prominent fixed point attracting 86% 
of all network states. The diagram shows the surprising 
feature that Pret{h) has its minimum at ft, = 1. This 
means that the main attractor, which is a fixed point, 
is not very stable under perturbations of one node. Fur- 
ther analysis shows that there are three other fixed points 
that differ from the main fixed point by the state of only 
one node. As suggested in [201, f^e cell is likely to be 
waiting for a new input when the network is at its fixed 
point. This might explain why this fixed point is rather 
sensitive to small perturbations. 
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FIG. 12: Pret (/i) (in percent) of the cell cycle regulation net- 
work of budding yeast 



V. CONCLUSIONS 

In this paper, we have evaluated the probability 
Pret{h) that a Boolean network returns to an attractor 
after perturbing h nodes, averaged over different initial 
states of the network. We found that Pret{h) can display 
a variety of different shapes, which yields insights into 
the state-space structure. If the response of the network 
to the perturbation of several nodes is independent for 
each perturbed node, Pret(h) decays exponentially with 
h for small h. Larger perturbations are of course not in- 
dependent, and if the perturbation leads the system to a 
random place in state space, Pret{h) shows a plateau for 
these perturbation sizes. When the size of the perturba- 
tion approaches the total number of nodes, Pretih) can 
shown various types of behaviour, including a decrease to 
zero or an increase to a value larger than the plateau. The 
reason is that the response to inverting the state of all 
nodes depends strongly on the network structure. We ob- 
tained analytical results for Random Boolean Networks 
in the limit of large A^. For critical networks, Pret{h) 
decreases rapidly to (almost) for large A^. For chaotic 
and frozen networks, Pret{h) remains nonzero for large 
A' and is not self-averaging, which means that the shape 
of Pret{h) differs widely between different networks. The 
ensemble average of the plateau value of Pret for chaotic 
networks is 2/3 for large A^. The fact that this plateau 
value is much higher than in critical networks means that 
chaotic networks are more robust than critical networks 
when perturbations affect a nonvanishing proportion of 
all nodes. The reason for this is that chaotic networks 
have a much smaller number of attractors than critical 
networks, and therefore the probability that a perturba- 
tion carries the network into the basin of attraction of 
a different attractor is larger in critical networks than 
in chaotic networks. Similarly interesting and surpris- 
ing is the result that a biological network that has been 
characterized as being very robust is pretty sensitive to 
perturbations of size /i = 1, while Pret is much larger for 
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larger perturbations. It remains to be seen if this feature VI. ACKNOWLEDGMENTS 

occurs also in other biological networks. 

To conclude, characterizing the dynamical stability of 
networks by a function Pret{h) gives far more informa- 
tion about the state space structure of the network than 

using a single number, such as the average "sensitivity" We thank C. Hamer for programming the yeast net- 

of nodes. work and T. Mihaljev for useful discussions. 
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